###Replication code - tables
#Electoral rewards and punishments for trade compensation####
#World Politics
#Minju KIM and Robert GULOTTY
#August 26, 2023

#####Preparation
library(stringr)
library(dplyr)
library(tidyr)
library(xtable)
library(tibble)
library(sandwich)
library(survival)
library(AER)
library(readr)
library(lubridate)
library(stargazer)
require(remotes)

###Packages needed to run the mediation analysis 
##Package "np": version 0.60.11
install_version("np", version = "0.60-11", repos = "http://cran.us.r-project.org")
library(np)
##Package "causalweight": version 0.1.0
install_version("causalweight", version = "0.1.0", repos = "http://cran.us.r-project.org")
library(causalweight)
##Package "glmnet": version 4.1.4 
library(glmnet)
##Package "mvtnorm": version 1.1.3
library(mvtnorm)
##Package "LARF": version 1.4
library(LARF)
##Package "hdm": version 0.3.1
library(hdm)
##Package "SuperLearner": version 2.0.28
library(SuperLearner)
##Package "xgboost": version 1.6.0.1
library(xgboost)
##Package "e1071": version 1.7.12 
library(e1071)
##Package "fastDummies": version 1.6.3
library(fastDummies)
##Package "grf": version 2.3.0
library(grf)
##Package "checkmate": version 2.1.0
library(checkmate)

select <- dplyr::select
mutate  <- dplyr::mutate
summarise <- dplyr::summarise
filter <- dplyr::filter

setwd("~/Dropbox/TAAChinaShock/World Politics/replication/")

#######################Tables in the manuscript
kpkg <- read.csv("czone_responsiveness.csv", as.is = T, na.strings = c("", NA))

kpkg <- kpkg %>% filter(year < 2007)

master <- read.csv(file = "master.csv", as.is = T)

###################################Table2: Is TAA responsive? By commuting zone
##Model1: Replicate Autor, Dorn, Hansen (ADH)
adh <- ivreg(
  d_trans_taaimp_pc ~ d_tradeusch_pw +
    l_shind_manuf_cbp + l_sh_routine33 +  l_task_outsource +
    l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn + t2 +
    reg_midatl + reg_encen + reg_wncen + reg_satl +
    reg_escen + reg_wscen + reg_mount +  reg_pacif |
    d_tradeotch_pw_lag + l_shind_manuf_cbp + l_sh_routine33 +  l_task_outsource +
    l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn + t2 +
    reg_midatl + reg_encen + reg_wncen + reg_satl +
    reg_escen + reg_wscen + reg_mount +  reg_pacif,
  data = kpkg,
  weights = timepwt48
)

m1 <- coeftest(adh, vcovCL, cluster = kpkg$statefip)


##Model2-Model3: Replicate Kim & Pelc (2021)
kp_cert_petit <- ivreg(
  d_cz_TAAdlrs_petit_cert_czpc ~
    d_tradeusch_pw + l_shind_manuf_cbp + l_sh_routine33 +
    l_task_outsource +
    l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn + t2 +
    reg_midatl + reg_encen + reg_wncen + reg_satl +
    reg_escen + reg_wscen + reg_mount +  reg_pacif |
    d_tradeotch_pw_lag + l_shind_manuf_cbp +
    l_sh_routine33 +  l_task_outsource +
    l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn + t2 +
    reg_midatl + reg_encen + reg_wncen + reg_satl +
    reg_escen + reg_wscen + reg_mount +  reg_pacif,
  data = kpkg,
  weights = timepwt48
)

m2 <- coeftest(kp_cert_petit, vcovCL, cluster = kpkg$statefip)


kp_cert_wrks <- ivreg(
  d_cz_TAAdlrs_wrks_cert_czpc ~
    d_tradeusch_pw + l_shind_manuf_cbp +
    l_sh_routine33 +  l_task_outsource +
    l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn + t2 +
    reg_midatl + reg_encen + reg_wncen + reg_satl +
    reg_escen + reg_wscen + reg_mount +  reg_pacif |
    d_tradeotch_pw_lag + l_shind_manuf_cbp + l_sh_routine33 +
    l_task_outsource +
    l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn + t2 +
    reg_midatl + reg_encen + reg_wncen + reg_satl +
    reg_escen + reg_wscen + reg_mount +  reg_pacif,
  data = kpkg,
  weights = timepwt48
)

m3 <- coeftest(kp_cert_wrks, vcovCL, cluster = kpkg$statefip)


##Model4-Model5: Kim and Gulotty (2023)
kg_cert_petit_n <- ivreg(
  d_cz_petit_cert ~
    d_tradeusch_pw + l_shind_manuf_cbp +
    l_sh_routine33 +  l_task_outsource +
    l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn + t2 +
    reg_midatl + reg_encen + reg_wncen + reg_satl +
    reg_escen + reg_wscen + reg_mount +  reg_pacif |
    d_tradeotch_pw_lag + l_shind_manuf_cbp +
    l_sh_routine33 + l_task_outsource +
    l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn + t2 +
    reg_midatl + reg_encen + reg_wncen + reg_satl +
    reg_escen + reg_wscen + reg_mount +  reg_pacif,
  data = kpkg,
  weights = timepwt48
)

m4 <- coeftest(kg_cert_petit_n, vcovCL, cluster = kpkg$statefip)


kg_cert_wrks_n <- ivreg(
  d_cz_wrks_cert ~
    d_tradeusch_pw + l_shind_manuf_cbp +
    l_sh_routine33 + l_task_outsource +
    l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn + t2 +
    reg_midatl + reg_encen + reg_wncen + reg_satl +
    reg_escen + reg_wscen + reg_mount +  reg_pacif |
    d_tradeotch_pw_lag + l_shind_manuf_cbp +
    l_sh_routine33 +  l_task_outsource +
    l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn + t2 +
    reg_midatl + reg_encen + reg_wncen + reg_satl +
    reg_escen + reg_wscen + reg_mount +  reg_pacif,
  data = kpkg,
  weights = timepwt48
)

m5 <- coeftest(kg_cert_wrks_n, vcovCL, cluster = kpkg$statefip)


stargazer(
  m1,
  m2,
  m3,
  m4,
  m5,
  type="html", out="~/Dropbox/TAAChinaShock/World Politics/replication/tables_replicated/table2.doc",
  title = "Is TAA responsive? By Commuting Zone",
  model.names = FALSE,
  notes = "N=722 CZs, excluding CZs in AK and HI",
  covariate.labels = c("$\\Delta$ Import Penetration"),
  dep.var.labels = c("Ten-Year Equivalent Change in TAA Distribution"),
  column.labels = c("ADH", "KP", "KG"),
  column.separate = c(1, 2, 2),
  add.lines = list(
    c("Observations", "1,444", "1,444", "1,444", "1,444", "1444"),
    c("Estimation", "2SLS", "2SLS", "2SLS", "2SLS", "2SLS"),
    c("2000 Ind/Occ Controls", "Yes", "Yes", "Yes", "Yes", "Yes"),
    c("2000 Demography Controls", "Yes", "Yes", "Yes", "Yes", "Yes"),
    c("Decade FE", "Yes", "Yes", "Yes", "Yes", "Yes")
  ),
  omit = c("l_", "t2", "reg_", "Constant"),
  omit.stat = c("rsq", "adj.rsq", "ser", "f"),
  keep.stat = c("n"),
  digits = 2,
  font.size = "small"
)



#########Table3: Heterogeneous Effect of TAA on Changes in Republican Vote Share
model1_int <- lm(d_rep_0800 ~ d_tradeusch_pw * pcert + t2,
                 weights = votetotal_2000,
                 data = master)

model2_int <-
  ivreg(
    d_rep_0800 ~ d_tradeusch_pw * pcert + t2 |
      d_tradeotch_pw_lag * pcert + t2,
    weights = votetotal_2000,
    data = master
  )

#Model 2: 1st stage F stat: 1236.4
first_stage <- lm(
  d_tradeusch_pw * pcert ~ t2 +
    d_tradeotch_pw_lag * pcert,
  weights = votetotal_2000,
  data = master
)
instrFtest <- waldtest(first_stage, . ~ . - d_tradeotch_pw_lag * pcert)
print(instrFtest)


model3_int <- ivreg(
  d_rep_0800 ~ d_tradeusch_pw * pcert + t2 +
    l_shind_manuf_cbp + l_sh_routine33 +  l_task_outsource |
    d_tradeotch_pw_lag * pcert + t2 +
    l_shind_manuf_cbp + l_sh_routine33 +  l_task_outsource,
  weights = votetotal_2000,
  data = master
)

#Model 3: 1st stage F stat: 645.56
first_stage <- lm(
  d_tradeusch_pw * pcert ~ t2 +
    l_shind_manuf_cbp + l_sh_routine33 +  l_task_outsource +
    d_tradeotch_pw_lag * pcert,
  weights = votetotal_2000,
  data = master
)
instrFtest <- waldtest(first_stage, . ~ . - d_tradeotch_pw_lag * pcert)
print(instrFtest)


model4_int <- ivreg(
  d_rep_0800 ~ d_tradeusch_pw * pcert + t2 +
    l_shind_manuf_cbp + l_sh_routine33 + l_task_outsource +
    l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn |
    d_tradeotch_pw_lag * pcert + t2 +
    l_shind_manuf_cbp + l_sh_routine33 +  l_task_outsource +
    l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn,
  weights = votetotal_2000,
  data = master
)

#Model 4: 1st stage F stat: 627.05
first_stage <- lm(
  d_tradeusch_pw * pcert ~ t2 +
    l_shind_manuf_cbp + l_sh_routine33 +  l_task_outsource++l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn +
    d_tradeotch_pw_lag * pcert,
  weights = votetotal_2000,
  data = master
)
instrFtest <- waldtest(first_stage, . ~ . - d_tradeotch_pw_lag * pcert)
print(instrFtest)

model5_int <- ivreg(
  d_rep_0800 ~ d_tradeusch_pw * pcert + t2 +
    l_shind_manuf_cbp + l_sh_routine33 +  l_task_outsource +
    l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn +
    region |
    d_tradeotch_pw_lag * pcert + t2 +
    l_shind_manuf_cbp + l_sh_routine33 +  l_task_outsource +
    l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn +
    region,
  weights = votetotal_2000,
  data = master
)

#Model 5: 1st stage F stat: 600.52
first_stage <- lm(
  d_tradeusch_pw * pcert ~ t2 +
    l_shind_manuf_cbp + l_sh_routine33 +  l_task_outsource++l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn +
    region +
    d_tradeotch_pw_lag * pcert,
  weights = votetotal_2000,
  data = master
)
instrFtest <- waldtest(first_stage, . ~ . - d_tradeotch_pw_lag * pcert)
print(instrFtest)



stargazer(
  model1_int,
  model2_int,
  model3_int,
  model4_int,
  model5_int,
  type="html", out="~/Dropbox/TAAChinaShock/World Politics/replication/tables_replicated/table3.doc",
  title = "Heterogeneous Effect of TAA on Changes in Republican Vote Share",
  model.names = FALSE,
  notes = "N=2,499 counties, excluding counties in AK and HI",
  covariate.labels = c(
    "$\\Delta$ Import Penetration",
    "Instrumented TAA Cert. Rate",
    "$\\Delta$ Import Penetratin* Instrumented TAA Cert. Rate"
  ),
  dep.var.labels = c("Changes in Republican Party Vote Share"),
  column.labels = c("OLS", "2SLS", "+Ind/Occ", "+Demography", "+Region"),
  add.lines = list(
    c("Estimation", "OLS", "2SLS", "2SLS", "2SLS", "2SLS"),
    c("2000 Ind/Occ Controls", "No", "No", "Yes", "Yes", "Yes"),
    c("2000 Demography Controls", "No", "No", "No", "Yes", "Yes"),
    c("Census Division Dummies", "No", "No", "No", "No", "Yes"),
    c("1st Stage F-Stat", "", "1236.4", "645.56", "627.05", "600.52")
  ),
  omit = c("l_", "t2", "region"),
  omit.stat = c("rsq", "adj.rsq", "ser", "f"),
  notes.append = FALSE,
  digits = 2,
  font.size = "small"
)



#################################Table4: Instrumented Causal Mediation Analysis
#####For using causalweight, only select variables that are needed in the dataset
master <- read.csv("master.csv")

master_med <- master%>%
  select(d_rep_0800, ipw_us2, actualcert, 
         ipw_others2, pcert, t2, l_shind_manuf_cbp, l_sh_routine33,
         l_task_outsource, l_sh_empl_f, l_sh_popedu_c, 
         l_sh_popfborn, Northeast, South, West, votetotal_2000)%>% 
  na.omit()


###Caution: The code below runs for a very long time.
#Note that the estimates are based on a simulation. 
set.seed(2948572)

mediation_500 <- medlateweight(
  y = master_med$d_rep_0800,
  d = master_med$ipw_us2,
  m = master_med$actualcert,
  zd = master_med$ipw_others2,
  zm = master_med$pcert,
  x = master_med$l_shind_manuf_cbp +
    master_med$t2 + master_med$l_sh_routine33  +
    master_med$l_task_outsource +
    master_med$l_sh_empl_f + master_med$l_sh_popedu_c +
    master_med$l_sh_popfborn +
    master_med$Northeast + master_med$South + master_med$West,
  boot = 500
)

mediation_500$results


stargazer(
  mediation_500$results,
  type="html",
  out="~/Dropbox/TAAChinaShock/World Politics/replication/tables_replicated/table4.doc",
  title = "Instrumented Causal Mediation Analysis",
  digits = 2,
  font.size = "small",
  style = "apsr"
)


#######################Tables in the appendix
######Table S1: Is TAA Responsive? (DV:Decade-average Change in TAA Distribution)
#TAA data
taazipcz_vote <- read.csv("taazipcz_vote.csv", as.is = T, 
                          na.strings = c("", "NA"))

workfile <- read.dta13("workfile_china.dta")

#Extract the distinct list of CZs from the China Shock dataset: 722 CZs 
czlist <- workfile%>%select(czone)
czlist <- distinct(czlist)

taa <- taazipcz_vote %>%
  filter(D_year > 1989 & D_year < 2008) %>%
  mutate(period2 = ifelse(D_year > 1999, 1, 0),
         certified_wrks = as.numeric(ifelse(certified == 1, Est..No..Workers, 0)))

taa_p1 <- taa %>% filter(period2 == 0)
taa_p2 <- taa %>% filter(period2 == 1)

bycz_p1 <- taa_p1 %>% group_by(czid_1990) %>%
  summarise(
    mcertified_p1 = sum(certified, na.rm = T) / 10,
    mwrks_p1 = sum(certified_wrks, na.rm = T) / 10
  )

bycz_p1_adj <- czlist %>%
  left_join(bycz_p1, by = c("czone" = "czid_1990")) %>%
  mutate_at(c(2:3), ~ replace(., is.na(.), 0))

bycz_p2 <- taa_p2 %>% group_by(czid_1990) %>%
  summarise(
    mcertified_p2 = sum(certified, na.rm = T) / 8,
    mwrks_p2 = sum(certified_wrks, na.rm = T) / 8
  )

bycz_p2_adj <- czlist %>%
  left_join(bycz_p2, by = c("czone" = "czid_1990")) %>%
  mutate_at(c(2:3), ~ replace(., is.na(.), 0))

bycz_p1p2 <- bycz_p1_adj %>%
  left_join(bycz_p2_adj) %>%
  mutate(d_mcertified = mcertified_p2 - mcertified_p1,
         d_mwrks = mwrks_p2 - mwrks_p1)

#China Shock data
workfile2 <-
  workfile %>% select(
    czone,
    statefip,
    city,
    yr,
    t2,
    d_tradeusch_pw,
    d_tradeotch_pw_lag,
    l_shind_manuf_cbp,
    l_sh_routine33,
    l_task_outsource,
    l_sh_empl_f,
    l_sh_popedu_c,
    l_sh_popfborn,
    l_sh_empl
  )

#Combine the TAA and China Shock data
combined <- bycz_p1p2 %>%
  left_join(workfile2, by = c("czone")) %>%
  arrange(czone)

###Let's run the regressions with the decade average
cz_cert_controls_final2 <-
  ivreg(
    d_mcertified ~ d_tradeusch_pw + l_shind_manuf_cbp + l_sh_routine33 +  l_task_outsource +
      l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn + t2 |
      d_tradeotch_pw_lag + l_shind_manuf_cbp + l_sh_routine33 +  l_task_outsource +
      l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn + t2,
    data = combined
  )

cz_cert_wrks_controls_final2 <-
  ivreg(
    d_mwrks ~ d_tradeusch_pw + l_shind_manuf_cbp + l_sh_routine33 +  l_task_outsource +
      l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn + t2 |
      d_tradeotch_pw_lag + l_shind_manuf_cbp + l_sh_routine33 +  l_task_outsource +
      l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn + t2,
    data = combined
  )


a <-
  coeftest(cz_cert_controls_final2, vcovCL, cluster = combined$statefip)
b <-
  coeftest(cz_cert_wrks_controls_final2, vcovCL, cluster = combined$statefip)

stargazer(
  a,
  b,
  type="html", out="~/Dropbox/TAAChinaShock/World Politics/replication/tables_replicated/tableS1.doc",
  title = "Is TAA responsive? By Commuting Zone",
  model.names = FALSE,
  notes = "N=722 CZs, excluding CZs in AK and HI",
  covariate.labels = c("$\\Delta$ Import Penetration"),
  dep.var.labels = c("Decade-average Change in TAA Distribution"),
  column.labels = c("Petitions", "Workers"),
  add.lines = list(
    c("Observations", "1,444", "1,444"),
    c("Estimation", "2SLS", "2SLS"),
    c("2000 Ind/Occ Controls", "Yes", "Yes"),
    c("2000 Demography Controls", "Yes", "Yes"),
    c("Decade FE", "Yes", "Yes")
  ),
  omit = c("l_", "t2", "Constant"),
  omit.stat = c("rsq", "adj.rsq", "ser", "f"),
  keep.stat = c("n"),
  digits = 2,
  font.size = "small"
)


######Table S2: Is TAA Responsive? (DV: Ten-year Change per person)
#First, calculate the ten-year change
taa2 <-
  taazipcz_vote %>% filter(D_year == 1990 |
                             D_year == 2000 | D_year == 2007) %>%
  mutate(certified_wrks = as.numeric(ifelse(certified == 1, Est..No..Workers, 0)))

bycz2 <- taa2 %>% group_by(czid_1990, D_year) %>%
  summarise(
    total = n(),
    numwrks = sum(as.numeric(Est..No..Workers), na.rm = T),
    certified = sum(certified),
    certified_numwrks = sum(certified_wrks, na.rm = T),
    leniency = certified / total * 100
  )

#Extract the distinct list of CZs from the China Shock dataset: 722 CZs (same as ADH, Kim and Pelc)
czlist <- workfile2 %>% select(czone)
czlist <- distinct(czlist)

#################################Rearrange the dataset at the level of cz
#The 1990 data
bycz_1990 <- bycz2 %>% filter(D_year == 1990) %>%
  select(czid_1990, total, numwrks, certified, certified_numwrks) %>%
  rename(
    total90 = total,
    numwrks90 = numwrks,
    certified90 = certified,
    certified_numwrks90 = certified_numwrks
  )

#Fill in CZs with zero TAA usage with zeros
bycz_1990_2 <- czlist %>%
  left_join(bycz_1990, by = c("czone" = "czid_1990")) %>%
  mutate_at(c(2:5), ~ replace(., is.na(.), 0))


#The 2000 data
bycz_2000 <- bycz2 %>% filter(D_year == 2000) %>%
  select(czid_1990, total, numwrks, certified, certified_numwrks) %>%
  rename(
    total00 = total,
    numwrks00 = numwrks,
    certified00 = certified,
    certified_numwrks00 = certified_numwrks
  )

#Fill in CZs with zero TAA usage with zeros
bycz_2000_2 <- czlist %>%
  left_join(bycz_2000, by = c("czone" = "czid_1990")) %>%
  mutate_at(c(2:5), ~ replace(., is.na(.), 0))


#The 2007 data
bycz_2007 <- bycz2 %>% filter(D_year == 2007) %>%
  select(czid_1990, total, numwrks, certified, certified_numwrks) %>%
  rename(
    total07 = total,
    numwrks07 = numwrks,
    certified07 = certified,
    certified_numwrks07 = certified_numwrks
  )

#Fill in CZs with zero TAA usage with zeros
bycz_2007_2 <- czlist %>%
  left_join(bycz_2007, by = c("czone" = "czid_1990")) %>%
  mutate_at(c(2:5), ~ replace(., is.na(.), 0))

#Combine the 1990, 2000, 2007 data: Calculate the difference in differences
bycz_d <-
  bycz_1990_2 %>% left_join(bycz_2000_2) %>% left_join(bycz_2007_2) %>%
  mutate(
    d_cert_0090 = certified00 - certified90,
    d_cert_0700 = certified07 - certified00,
    d_cert_wrks_0090 = certified_numwrks00 - certified_numwrks90,
    d_cert_wrks_0700 = certified_numwrks07 - certified_numwrks00,
    did_cert = d_cert_0700 - d_cert_0090,
    did_cert_wrks = d_cert_wrks_0700 - d_cert_wrks_0090,
    did_cert_adj = d_cert_0700 - d_cert_0090 * 10 / 7,
    did_cert_wrks_adj = d_cert_wrks_0700 - d_cert_wrks_0090 * 10 /
      7
  )

bycz_d_0090 <- bycz_d %>%
  select(czone, d_cert_0090, d_cert_wrks_0090) %>%
  mutate(t2 = 0) %>%
  rename(d_cert = d_cert_0090, d_cert_wrks = d_cert_wrks_0090)

bycz_d_0700 <- bycz_d %>%
  select(czone, d_cert_0700, d_cert_wrks_0700) %>%
  mutate(t2 = 1) %>%
  rename(d_cert = d_cert_0700, d_cert_wrks = d_cert_wrks_0700)

combined <- rbind(bycz_d_0090, bycz_d_0700)

combined2 <- combined %>%
  left_join(workfile2, by = c("czone", "t2")) %>%
  arrange(czone)

czpop <- read.csv("cz2000_DOA.csv")

czpop <-
  unique(czpop %>% select(czid_1990, cz_ppl_2000_doa)) %>% filter(complete.cases(.))
czpop$cz_ppl_2000_doa <-
  as.numeric(gsub(",", "", czpop$cz_ppl_2000_doa))
#The data is basedo n ppl in cz_2000, so let's combine in case cz in 1990 was split into two
czpop_adj <- czpop %>% group_by(czid_1990) %>%
  summarise(cz_ppl_1990 = sum(cz_ppl_2000_doa))

combined2 <-
  combined2 %>% left_join(czpop_adj, by = c("czone" = "czid_1990"))

combined2 <- combined2 %>%
  mutate(
    d_cert_pw = d_cert / cz_ppl_1990,
    d_cert_wrks_pw = d_cert_wrks / cz_ppl_1990,
    d_cert_pw_empl = d_cert / (cz_ppl_1990 * l_sh_empl / 100),
    d_cert_wrks_pw_empl = d_cert_wrks / (cz_ppl_1990 * l_sh_empl /
                                           100)
  )


###Let's run the regressions with the decade average per person
cz_cert_controls_final3 <-
  ivreg(
    d_cert_pw ~ d_tradeusch_pw + l_shind_manuf_cbp + l_sh_routine33 +  l_task_outsource +
      l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn + t2 |
      d_tradeotch_pw_lag + l_shind_manuf_cbp + l_sh_routine33 +  l_task_outsource +
      l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn + t2,
    data = combined2
  )

cz_cert_wrks_controls_final3 <-
  ivreg(
    d_cert_wrks_pw ~ d_tradeusch_pw + l_shind_manuf_cbp + l_sh_routine33 +  l_task_outsource +
      l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn + t2
    |
      d_tradeotch_pw_lag + l_shind_manuf_cbp + l_sh_routine33 +  l_task_outsource +
      l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn + t2,
    data = combined2
  )


c <-
  coeftest(cz_cert_controls_final3, vcovCL, cluster = combined2$statefip)

d <-
  coeftest(cz_cert_wrks_controls_final3, vcovCL, cluster = combined2$statefip)


stargazer(
  c,
  d,
  type="html", out="~/Dropbox/TAAChinaShock/World Politics/replication/tables_replicated/tableS2.doc",
  title = "Is TAA responsive? By Commuting Zone",
  model.names = FALSE,
  notes = "N=722 CZs, excluding CZs in AK and HI",
  covariate.labels = c("$\\Delta$ Import Penetration"),
  dep.var.labels = c("Ten-Year Change per person"),
  column.labels = c("Petitions", "Workers"),
  add.lines = list(
    c("Observations", "1,444", "1,444"),
    c("Estimation", "2SLS", "2SLS"),
    c("2000 Ind/Occ Controls", "Yes", "Yes"),
    c("2000 Demography Controls", "Yes", "Yes"),
    c("Decade FE", "Yes", "Yes")
  ),
  omit = c("l_", "t2", "Constant"),
  omit.stat = c("rsq", "adj.rsq", "ser", "f"),
  keep.stat = c("n"),
  digits = 6,
  font.size = "small"
)


######Table S3: Is TAA Responsive? (DV: Ten-year Change per worker)
cz_cert_controls_final4 <-
  ivreg(
    d_cert_pw_empl ~ d_tradeusch_pw + l_shind_manuf_cbp + l_sh_routine33 +  l_task_outsource +
      l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn + t2 |
      d_tradeotch_pw_lag + l_shind_manuf_cbp + l_sh_routine33 +  l_task_outsource +
      l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn + t2,
    data = combined2
  )

cz_cert_wrks_controls_final4 <-
  ivreg(
    d_cert_wrks_pw_empl ~ d_tradeusch_pw + l_shind_manuf_cbp + l_sh_routine33 +  l_task_outsource +
      l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn + t2
    |
      d_tradeotch_pw_lag + l_shind_manuf_cbp + l_sh_routine33 +  l_task_outsource +
      l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn + t2,
    data = combined2
  )


e <-
  coeftest(cz_cert_controls_final4, vcovCL, cluster = combined2$statefip)

f <-
  coeftest(cz_cert_wrks_controls_final4, vcovCL, cluster = combined2$statefip)

stargazer(
  e,
  f,
  type="html", out="~/Dropbox/TAAChinaShock/World Politics/replication/tables_replicated/tableS3.doc",
  title = "Is TAA responsive? By Commuting Zone",
  model.names = FALSE,
  notes = "N=722 CZs, excluding CZs in AK and HI",
  covariate.labels = c("$\\Delta$ Import Penetration"),
  dep.var.labels = c("Ten-Year Change in TAA Distribution, per worker"),
  column.labels = c("Petitions", "Petitions", "Workers", "Workers"),
  add.lines = list(
    c("Observations", "1,444", "1,444"),
    c("Estimation", "2SLS", "2SLS"),
    c("2000 Ind/Occ Controls", "Yes", "Yes"),
    c("2000 Demography Controls", "Yes", "Yes"),
    c("Decade FE", "Yes", "Yes")
  ),
  omit = c("l_", "t2", "Constant"),
  omit.stat = c("rsq", "adj.rsq", "ser", "f"),
  keep.stat = c("n"),
  digits = 5,
  font.size = "small"
)


######Table S4: Summary statistics
master_ss <- master %>% select(
  d_rep_0800,
  d_tradeusch_pw,
  d_tradeotch_pw_lag,
  ipw_us2,
  ipw_others2,
  l_shind_manuf_cbp,
  l_sh_routine33,
  l_task_outsource,
  l_sh_empl_f,
  l_sh_popedu_c,
  l_sh_popfborn,
  Midwest,
  Northeast,
  South,
  West,
  pcert,
  actualcert,
  totalpetitions
)

stargazer(
  master_ss,
  type="html", out="~/Dropbox/TAAChinaShock/World Politics/replication/tables_replicated/tableS4.doc",
  nobs = FALSE,
  mean.sd = TRUE,
  median = TRUE,
  title =  "Summary Statistics (N=4,998)",
  covariate.labels = c(
    "$\\Delta$Rep 2008-2000",
    "$\\Delta$US-China trade pw",
    "$\\Delta$Lagged-Other trade pw",
    "Binary US-China Trade",
    "Binary Other Trade",
    "Lagged Share Manufacturing",
    "Lagged Share Routine",
    "Lagged Offshorability",
    "Lagged Share Female Employed",
    "Lagged Share Collage Population",
    "Lagged Share Foreign Born",
    "Midwest",
    "Northeast",
    "South",
    "West",
    "Predicted TAA Certification",
    "Actual TAA Certification",
    "Total TAA Certification"
  ),
  digits = 2,
  font.size = "footnotesize"
)


######Table S5: Statistics of TAA Certifying Officers
#####TAA from 2005 to 2008
taazipcz_cut <- taazipcz_vote %>%
  select(D_year, untenured, certified, Det.Date, Cert.Officer, czid_1990) %>%
  filter(Det.Date > as.Date("2004-11-02") &
           Det.Date < as.Date("2008-11-04"))

sum(is.na(taazipcz_cut$Cert.Officer))
###46 NAs out of 9636 observations (0.005%)

##eliminate the NAs (left with 9035 observations)
taazipcz_cut <- taazipcz_cut %>%
  filter(!is.na(certified) & !is.na(Cert.Officer) & !is.na(czid_1990))

###Certifying Officer Statistics
costat <- taazipcz_cut %>% group_by(Cert.Officer) %>%
  summarise(
    total = n(),
    leniency = mean(certified),
    firstcase = min(Det.Date),
    lastcase = max(Det.Date)
  )

costat_df <- as.tibble(costat)
costat_df <- xtable(costat_df, 
       caption = "Statistics of TAA Certifying Officers", 
       digits =2)

print(costat_df, 
      type="html", 
      file="~/Dropbox/TAAChinaShock/World Politics/replication/tables_replicated/tableS5.doc")

######Table S6: Random Assignment of TAA petitions to Certifying Officers
master2 <- master %>%
  filter(!is.na(pcert) & !is.na(region)) %>%
  filter(t2 == 0) %>%
  select(
    czid_1990,
    l_shind_manuf_cbp,
    l_sh_routine33,
    l_task_outsource,
    l_sh_empl_f,
    l_sh_popedu_c,
    l_sh_popfborn,
    Midwest,
    Northeast,
    South,
    West
  )

#Insert the predicted certification rate of each certifying officer
taazipcz_cut2 <- taazipcz_cut %>%
  left_join(master2) %>%
  unique() %>%
  mutate(
    Cert_leniency = case_when(
      Cert.Officer == "Kushner" ~ 0.6773050,
      Cert.Officer == "Poole" ~ 0.6416868,
      Cert.Officer == "Church" ~ 0.5808246
    )
  )

b <- lm(l_sh_empl_f ~ Cert.Officer, data = taazipcz_cut2)
c <- lm(l_sh_popedu_c ~ Cert.Officer, data = taazipcz_cut2)
d <- lm(l_shind_manuf_cbp ~ Cert.Officer, data = taazipcz_cut2)
e <- lm(l_sh_popfborn ~ Cert.Officer, data = taazipcz_cut2)
f <- lm(l_sh_routine33 ~ Cert.Officer, data = taazipcz_cut2)
g <- lm(l_task_outsource ~ Cert.Officer, data = taazipcz_cut2)
h <- lm(Midwest ~ Cert.Officer, data = taazipcz_cut2)
i <- lm(Northeast ~ Cert.Officer, data = taazipcz_cut2)
k <- lm(West ~ Cert.Officer, data = taazipcz_cut2)

ag <- rbind(
  tidy(b)[2:3,],
  tidy(c)[2:3,],
  tidy(d)[2:3,],
  tidy(e)[2:3,],
  tidy(f)[2:3,],
  tidy(g)[2:3,],
  tidy(h)[2:3,],
  tidy(i)[2:3,],
  tidy(k)[2:3,]
) %>%
  select(-statistic) %>%
  mutate(
    variable =  c(
      "Lagged Share of Female Employed",
      "Lagged Share of Female Employed",
      "Lagged Share of College Population",
      "Lagged Share of College Population",
      "Lagged  Share of Manufacturing",
      "Lagged  Share of Manufacturing",
      "Lagged Share of Foreign Born",
      "Lagged Share of Foreign Born",
      "Lagged Share of Routineness",
      "Lagged Share of Routineness",
      "Lagged Offshorability",
      "Lagged Offshorability",
      "Midwest",
      "Midwest",
      "Northwest",
      "Northwest",
      "West",
      "West"
    ),
    cert.officer = case_when(
      term == "Cert.OfficerKushner" ~ "B",
      term == "Cert.OfficerPoole" ~ "C"
    )
  )
ag <- ag %>%
  select(variable, cert.officer, estimate, std.error, p.value)

breg <- xtable(ag, caption = "Random Assignment of CZ-level Predicted TAA Certification", 
         digits =2)

print(breg, type="html",
      file="~/Dropbox/TAAChinaShock/World Politics/replication/tables_replicated/tableS6.doc")

######Table S7:The Consternation Effect in TAA Certification
master$ipw_us2 <- as.factor(master$ipw_us2)
levels(master$ipw_us2) <- c("Low", "High")
master$ipw_others2 <- as.factor(master$ipw_others2)
levels(master$ipw_others2) <- c("Low", "High")

###Regression model with the instrument (with predicted TAA certification)
model_int_b <- ivreg(
  d_rep_0800 ~ ipw_us2 * pcert + t2 +
    l_shind_manuf_cbp + l_sh_routine33 +  l_task_outsource +
    l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn +
    region |
    ipw_others2 * pcert + t2 +
    l_shind_manuf_cbp + l_sh_routine33 +  l_task_outsource +
    l_sh_empl_f + l_sh_popedu_c + l_sh_popfborn +
    region,
  weights = votetotal_2000,
  data = master
)

stargazer(
  model_int_b,
  type="html",
  out="~/Dropbox/TAAChinaShock/World Politics/replication/tables_replicated/tableS7.doc",
  title = "The Consternation Effect in TAA Certification",
  model.names = FALSE,
  notes = "N=2,499 counties, excluding counties in AK and HI",
  covariate.labels = c(
    "High Economic Shock",
    "Instrumented Certification",
    "High Economic Shock*Instrumented Certification"
  ),
  dep.var.labels = c("Changes in Republican Party Vote Share"),
  add.lines = list(
    c("2000 Ind/Occ Controls", "Yes"),
    c("2000 Demography Controls", "Yes"),
    c("Census Division Dummies", "Yes")
  ),
  omit = c("l_", "t2", "region"),
  omit.stat = c("rsq", "adj.rsq", "ser", "f"),
  notes.append = FALSE,
  digits = 3,
  font.size = "normalsize"
)


######Table S8: Testing the Conditional Independence Assumption: Relationship between Chinese Exports to High-income Markets (DZ) and Instrumented TAA Certification (MZ)
master <- read.csv("master.csv")

dz_mz <-
  lm(
    ipw_others2 ~ pcert + l_sh_empl_f + l_sh_popedu_c + l_shind_manuf_cbp +
      l_sh_popfborn + l_sh_routine33 + l_task_outsource +
      Northeast + South + West,
    data = master
  )


stargazer(
  dz_mz,
  type="html",
  out="~/Dropbox/TAAChinaShock/World Politics/replication/tables_replicated/tableS8.doc",
  title = "Testing the Conditional Independence Assumption: Relationship between Chinese Exports to High-income Markets (DZ) and Instrumented TAA Certification (MZ)",
  covariate.labels = c(
    "Instrumented Certification Rate (MZ)",
    "Lagged Share of Female Employed",
    "Lagged Share of College Population",
    "Lagged  Share of Manufacturing",
    "Lagged Share of Foreign Born",
    "Lagged Share of Routineness",
    "Lagged Offshorability",
    "Northeest",
    "South",
    "West"
  ),
  dep.var.labels = c("Chinese Exports to High-income Markets (DZ)"),
  omit.stat = c("adj.rsq"),
  omit.table.layout = "n",
  digits = 3,
  font.size = "footnotesize"
)
